7 TYとの近接と被攻撃との関連

ここでは、交尾期においてTYやITとの近接がオスからの被攻撃頻度に及ぼす影響を調べていきます。

7.1 TY/ITとの近接と被攻撃頻度の関連

まずは、休息集団内にTYとITがいるか否かによって、メスの被攻撃頻度が変わるかを検討します。

7.1.1 非発情メス

非発情メスについての分析を行います。

7.1.1.1 データの加工

分析のためのデータを作成します。データの独立性を担保するため、3分ごとのデータを用います。そのほかは以下の通りです。

• 地上にいるデータのみを用いる
  
• 連続変数は標準化

## 地上のデータ、交尾期のデータを抽出  
focal_agg_prox_an <- focal_combined_prox_final %>% 
  filter(rs2 == 0) %>% 
  filter(T_G == "G") %>% 
  filter(time %% 3 == 0) %>% 
  filter(!str_detect(study_period, "nm")) %>% 
  mutate(focalID = str_c(study_period, "_", no_focal)) %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         x3m_female_std = standardize(x3m_female)) 

focal_agg_prox_an

7.1.1.2 モデリング

それでは、モデリングを行います。調査期間は多重共線性が大きくなるため除外しました。

• 応答変数: 直後1分間の攻撃の有無(agg_focal)
  
• 説明変数: 3m以内のメス数(x3m_female)、3m以内のオスの有無(TY_3m, IT_3m, KR_3m)、調査日の確認メス数(female_std)、調査日の発情メス数(est_std)、追跡中のTYとITの有無(TY_presence、IT_presence)
  
• ランダム切片: メスID(subject)、追跡セッションID(focalID)
  
• 分布: ベルヌーイ分布

## 全期間  
m_agg_prox_an <- brm(agg_focal ~ x3m_female_std + TY_3m + IT_3m + KR_3m +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_prox_an,
                       file = "model/m_agg_prox_an.rds")

## 2020年交尾期以前  
focal_agg_prox_an_bf21 <- focal_agg_prox_an %>% 
  filter(study_period != "m21") %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         x3m_female_std = standardize(x3m_female)) 

m_agg_prox_an_bf21 <- brm(agg_focal ~ x3m_female_std + TY_3m + IT_3m + KR_3m +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_prox_an_bf21,
                       file = "model/m_agg_prox_an_bf21.rds")

7.1.1.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_agg_prox_an <- dh_check_brms(m_agg_prox_an, quantreg = TRUE)

## 21年より前
dh_agg_prox_an_bf21 <- dh_check_brms(m_agg_prox_an_bf21, quantreg = TRUE)

多重共線性の問題もありません。

## 全期間   
check_collinearity(m_agg_prox_an)

## 21年より前
check_collinearity(m_agg_prox_an_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_agg_prox_an, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_agg_prox_an_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_agg_prox_an)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_agg_prox_an_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

最後に、残差の時系列相関を調べます。まずは全期間についてです。有意な時系列相関はないようです。

## 全時間ポイントの作成  
focal_agg_prox_an_full <-  focal_agg_prox_an %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_prox_an$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_prox_an_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_prox_an <- focal_agg_prox_an_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_prox_an_full2 <- focal_agg_prox_an_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_prox_an) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_prox_an_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_prox_an_full2$resid, type = "Ljung-Box", lag = 10)

## 
##  Box-Ljung test
## 
## data:  focal_agg_prox_an_full2$resid
## X-squared = 8.3965, df = 10, p-value = 0.5902

2021年交尾期を除いた場合も同じです。

## 全時間ポイントの作成  
focal_agg_prox_an_bf21_full <-  focal_agg_prox_an_bf21 %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_prox_an_bf21$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_prox_an_bf21_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_prox_an_bf21 <- focal_agg_prox_an_bf21_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_prox_an_bf21_full2 <- focal_agg_prox_an_bf21_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_prox_an_bf21) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_prox_an_bf21_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_prox_an_bf21_full2$resid, type = "Ljung-Box", lag = 10)

## 
##  Box-Ljung test
## 
## data:  focal_agg_prox_an_bf21_full2$resid
## X-squared = 12.85, df = 10, p-value = 0.2322

7.1.1.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。

続いて、2021年のデータを除いたモデルについて確認する。

7.1.2 発情メス

発情メスについての分析を行います。

7.1.2.1 データの加工

分析のためのデータを作成します。データの独立性を担保するため、3分ごとのデータを用います。そのほかは以下の通りです。

• 地上にいるデータのみを用いる
  
• 連続変数は標準化

## 地上のデータ、交尾期のデータを抽出  
focal_agg_prox_es <- focal_combined_prox_final %>% 
  filter(rs2 == 1) %>% 
  filter(T_G == "G") %>% 
  filter(time %% 3 == 0) %>% 
  filter(!str_detect(study_period, "nm")) %>% 
  mutate(focalID = str_c(study_period, "_", no_focal)) %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         x3m_female_std = standardize(x3m_female)) 

focal_agg_prox_es

7.1.2.2 モデリング

それでは、モデリングを行います。調査期間は多重共線性が大きくなるため除外しました。

• 応答変数: 直後1分間の攻撃の有無(agg_focal)
  
• 説明変数: 3m以内のメス数(x3m_female)、3m以内のオスの有無(TY_3m, IT_3m, KR_3m)、調査日の確認メス数(female_std)、調査日の発情メス数(est_std)、追跡中のTYとITの有無(TY_presence、IT_presence)
  
• ランダム切片: メスID(subject)、追跡セッションID(focalID)
  
• 分布: ベルヌーイ分布(21年を除いたものはベータ二項分布)

## 全期間  
m_agg_prox_es <- brm(agg_focal ~ x3m_female_std + TY_3m + IT_3m + KR_3m +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_prox_es,
                       file = "model/m_agg_prox_es.rds")

## 2020年交尾期以前  
focal_agg_prox_es_bf21 <- focal_agg_prox_es %>% 
  filter(study_period != "m21") %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         x3m_female_std = standardize(x3m_female)) %>% 
  mutate(rowid = 1:n())

## divergent transitionが少し出る  
m_agg_prox_es_bf21 <- brm(agg_focal|trials(1) ~ x3m_female_std + TY_3m + IT_3m + KR_3m +
                            female_std + est_std + TY_presence + IT_presence + 
                            (1|focalID) + (1|subject),
                          family = beta_binomial(),
                          iter = 11000, warmup = 1000, seed = 1106,
                          prior = c(prior(student_t(4,0,15), class = "b"),
                                    prior(student_t(4,0,15), class = "Intercept"),
                                    prior(student_t(4,0,10), class = "sd")),
                          control=list(adapt_delta = 0.99999, max_treedepth = 20),
                          backend = "cmdstanr",
                          data = focal_agg_prox_es_bf21,
                          file = "model/m_agg_prox_es_bf21.rds")

7.1.2.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。21年を除いた場合に少し悪そう。

## 全期間    
dh_agg_prox_es <- dh_check_brms(m_agg_prox_es, quantreg = TRUE)

## 21年より前
dh_agg_prox_es_bf21 <- dh_check_brms(m_agg_prox_es_bf21, quantreg = TRUE)

多重共線性の問題もありません。

## 全期間   
check_collinearity(m_agg_prox_es)

## 21年より前
check_collinearity(m_agg_prox_es_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_agg_prox_es, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_agg_prox_es_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_agg_prox_es)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_agg_prox_es_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

最後に、残差の時系列相関を調べます。まずは全期間についてです。有意な時系列相関はないようです。

## 全時間ポイントの作成  
focal_agg_prox_es_full <-  focal_agg_prox_es %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_prox_es$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_prox_es_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_prox_es <- focal_agg_prox_es_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_prox_es_full2 <- focal_agg_prox_es_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_prox_es) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_prox_es_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_prox_es_full2$resid, type = "Ljung-Box", lag = 10)

## 
##  Box-Ljung test
## 
## data:  focal_agg_prox_es_full2$resid
## X-squared = 3.1539, df = 10, p-value = 0.9776

2021年交尾期を除いた場合も同じです。

## 全時間ポイントの作成  
focal_agg_prox_es_bf21_full <-  focal_agg_prox_es_bf21 %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_prox_es_bf21$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_prox_es_bf21_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_prox_es_bf21 <- focal_agg_prox_es_bf21_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_prox_es_bf21_full2 <- focal_agg_prox_es_bf21_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_prox_es_bf21) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_prox_es_bf21_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_prox_es_bf21_full2$resid, type = "Ljung-Box", lag = 10)

## 
##  Box-Ljung test
## 
## data:  focal_agg_prox_es_bf21_full2$resid
## X-squared = 14.496, df = 10, p-value = 0.1515

7.1.2.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。

続いて、2021年のデータを除いたモデルについて確認する。

7.2 休息集団内のオスの有無と被攻撃頻度

7.2.1 非発情メス

非発情メスについての分析を行います。

7.2.1.1 データの加工

分析のためのデータを作成します。データの独立性を担保するため、3分ごとのデータを用います。そのほかは以下の通りです。

• 地上にいるデータのみを用いる
  
• 連続変数は標準化

## 地上のデータ、交尾期のデータを抽出  
focal_agg_RG_an <- focal_combined_prox_final %>% 
  filter(rs2 == 0) %>% 
  filter(RG == 1 & T_G == "G") %>% 
  filter(time %% 3 == 0) %>% 
  filter(!str_detect(study_period, "nm")) %>% 
  mutate(focalID = str_c(study_period, "_", no_focal)) %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         RG_female_std = standardize(RG_female)) 

7.2.1.2 モデリング

それでは、モデリングを行います。調査期間は多重共線性が大きくなるため除外しました。

• 応答変数: 直後1分間の攻撃の有無(agg_focal)
  
• 説明変数: 3m以内のメス数(x3m_female)、3m以内のオスの有無(TY_3m, IT_3m, KR_3m)、調査日の確認メス数(female_std)、調査日の発情メス数(est_std)、追跡中のTYとITの有無(TY_presence、IT_presence)
  
• ランダム切片: メスID(subject)、追跡セッションID(focalID)
  
• 分布: ベルヌーイ分布

## 全期間  
m_agg_RG_an <- brm(agg_focal ~ RG_female_std + RG_TY + RG_IT + RG_KR +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_RG_an,
                       file = "model/m_agg_RG_an.rds")

## 2020年交尾期以前  
focal_agg_RG_an_bf21 <- focal_agg_RG_an %>% 
  filter(study_period != "m21") %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         RG_female_std = standardize(RG_female)) 

m_agg_RG_an_bf21 <- brm(agg_focal ~ RG_female_std + RG_TY + RG_IT + RG_KR +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_RG_an_bf21,
                       file = "model/m_agg_RG_an_bf21.rds")

7.2.1.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_agg_RG_an <- dh_check_brms(m_agg_RG_an, quantreg = TRUE)

## 21年より前
dh_agg_RG_an_bf21 <- dh_check_brms(m_agg_RG_an_bf21, quantreg = TRUE)

多重共線性の問題もありません。

## 全期間   
check_collinearity(m_agg_RG_an)

## 21年より前
check_collinearity(m_agg_RG_an_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_agg_RG_an, ndraws = 100)+
  theme_bw() +
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_agg_RG_an_bf21, ndraws = 100)+
  theme_bw() +
  theme(aspect.ratio = 0.9)

Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_agg_RG_an)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_agg_RG_an_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

最後に、残差の時系列相関を調べます。まずは全期間についてです。有意な時系列相関はないようです。

## 全時間ポイントの作成  
focal_agg_RG_an_full <-  focal_agg_RG_an %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_RG_an$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_RG_an_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_RG_an <- focal_agg_RG_an_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_RG_an_full2 <- focal_agg_RG_an_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_RG_an) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_RG_an_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_RG_an_full2$resid, type = "Ljung-Box", lag = 10)

## 
##  Box-Ljung test
## 
## data:  focal_agg_RG_an_full2$resid
## X-squared = 15.81, df = 10, p-value = 0.1052

2021年交尾期を除いた場合も同じです。

## 全時間ポイントの作成  
focal_agg_RG_an_bf21_full <-  focal_agg_RG_an_bf21 %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_RG_an_bf21$scaledResiduals) 

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_RG_an_bf21_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

## フォーカル間に入れる空白行  
time_blank_agg_RG_an_bf21 <- focal_agg_RG_an_bf21_full %>% 
  group_by(focalID, date, subject) %>% 
  summarise(time = max(time) + 1) %>% 
  group_by(focalID, date, subject) %>% 
  complete(time = seq(time, time + 30, by = 3)) %>% 
  ungroup() %>% 
  mutate(resid = NA)

focal_agg_RG_an_bf21_full2 <- focal_agg_RG_an_bf21_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0) %>% 
  full_join(time_blank_agg_RG_an_bf21) %>% 
  arrange(focalID, time)

## 自己相関の図示  
acf(focal_agg_RG_an_bf21_full2$resid, lag.max = 10, na.action = na.pass)

## Box-Ljung test
Box.test(focal_agg_RG_an_bf21_full2$resid, type = "Ljung-Box", lag = 10)  

## 
##  Box-Ljung test
## 
## data:  focal_agg_RG_an_bf21_full2$resid
## X-squared = 29.385, df = 10, p-value = 0.001079

7.2.1.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。

続いて、2021年のデータを除いたモデルについて確認する。

7.2.2 発情メス

発情メスについての分析を行います。

7.2.2.1 データの加工

分析のためのデータを作成します。データの独立性を担保するため、3分ごとのデータを用います。そのほかは以下の通りです。

• 地上にいるデータのみを用いる
  
• 連続変数は標準化

## 地上のデータ、交尾期のデータを抽出  
focal_agg_RG_es <- focal_combined_prox_final %>% 
  filter(rs2 == 1) %>% 
  filter(RG == 1 & T_G == "G") %>% 
  filter(time %% 3 == 0) %>% 
  filter(!str_detect(study_period, "nm")) %>% 
  mutate(focalID = str_c(study_period, "_", no_focal)) %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         RG_female_std = standardize(RG_female)) 

7.2.2.2 モデリング

それでは、モデリングを行います。調査期間は多重共線性が大きくなるため除外しました。

• 応答変数: 直後1分間の攻撃の有無(agg_focal)
  
• 説明変数: 3m以内のメス数(x3m_female)、3m以内のオスの有無(TY_3m, IT_3m, KR_3m)、調査日の確認メス数(female_std)、調査日の発情メス数(est_std)、追跡中のTYとITの有無(TY_presence、IT_presence)
  
• ランダム切片: メスID(subject)、追跡セッションID(focalID)
  
• 分布: ベルヌーイ分布(21年を除いたものはベータ二項分布)

## 全期間  
m_agg_RG_es <- brm(agg_focal ~ RG_female_std + RG_TY + RG_IT + RG_KR +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_RG_es,
                       file = "model/m_agg_RG_es.rds")

## 2020年交尾期以前  
focal_agg_RG_es_bf21 <- focal_agg_RG_es %>% 
  filter(study_period != "m21") %>% 
  mutate(CSI_TY_std = standardize(CSI_TY),
         CSI_IT_std = standardize(CSI_IT),
         ntm_std = standardize(no_ntm),
         est_std = standardize(no_est),
         female_std = standardize(no_female),
         RG_female_std = standardize(RG_female)) %>% 
  mutate(rowid = 1:n())

m_agg_RG_es_bf21 <- brm(agg_focal|trials(1) ~ RG_female_std + RG_TY + RG_IT + RG_KR +
                         female_std + est_std + TY_presence + IT_presence + 
                         (1|focalID) + (1|subject),
                       family = binomial(),
                       iter = 11000, warmup = 1000, seed = 123,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.99999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = focal_agg_RG_es_bf21,
                       file = "model/m_agg_RG_es_bf21.rds")

7.2.2.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_agg_RG_es <- dh_check_brms(m_agg_RG_es, quantreg = TRUE)

## 21年より前
dh_agg_RG_es_bf21 <- dh_check_brms(m_agg_RG_es_bf21, quantreg = TRUE)

多重共線性の問題もありません。

## 全期間   
check_collinearity(m_agg_RG_es)

## 21年より前
check_collinearity(m_agg_RG_es_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_agg_RG_es, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_agg_RG_es_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_agg_RG_es)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_agg_RG_es_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

最後に、残差の時系列相関を調べます。まずは全期間についてです。有意な時系列相関はないようです。

## 全時間ポイントの作成  
focal_agg_RG_es_full <-  focal_agg_RG_es%>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_RG_es$scaledResiduals) %>% 
  group_by(focalID, subject, date) %>% 
  complete(time = full_seq(time, 1)) %>% 
  ungroup()

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_RG_es_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

focal_agg_RG_es_full2 <- focal_agg_RG_es_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0)
  
testTemporalAutocorrelation(dh_agg_RG_es, time = drop_na(focal_agg_RG_es_full2, 
                                                           agg_focal) %>% .$actualtime)

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.013, p-value = 0.7737
## alternative hypothesis: true autocorrelation is not 0

2021年交尾期を除いた場合も同じです。

## 全時間ポイントの作成  
focal_agg_RG_es_bf21_full <-  focal_agg_RG_es_bf21 %>% 
  drop_na(agg_focal) %>% 
  mutate(resid = dh_agg_RG_es_bf21$scaledResiduals) %>% 
  group_by(focalID, subject, date) %>% 
  complete(time = full_seq(time, 1)) %>% 
  ungroup()

## 日ごとに個体追跡セッション番号を作成  
date_focal_list <- focal_agg_RG_es_bf21_full %>% 
  distinct(date, focalID) %>% 
  group_by(date) %>% 
  mutate(dateID = 1:n()) %>% 
  ungroup() 

focal_agg_RG_es_bf21_full2 <- focal_agg_RG_es_bf21_full %>% 
  left_join(date_focal_list) %>% 
  ## 時刻を作成
  ## 2018年はかぶらないように作成  
  group_by(focalID) %>% 
  mutate(start_time = case_when(
    study_period == "m18" ~ make_datetime(year(date), month(date), mday(date),
                                          hour = unique(dateID*3), min = 0, sec = 0),
    .default = start_time
  )) %>% 
  ungroup() %>% 
  mutate(actualtime = start_time + time) %>% 
  group_by(focalID) %>% 
  complete(time = full_seq(time, 1)) %>% 
  filter(time %% 3 == 0)
  
testTemporalAutocorrelation(dh_agg_RG_es_bf21, time = drop_na(focal_agg_RG_es_bf21_full2, 
                                                           agg_focal) %>% .$actualtime)

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 2.1155, p-value = 0.05199
## alternative hypothesis: true autocorrelation is not 0

7.2.2.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。

続いて、2021年のデータを除いたモデルについて確認する。

References

Gabry, J., & Mahr, T. (2022). Bayesplot: Plotting for bayesian models. https://mc-stan.org/bayesplot

Hartig, F. (2022). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. https://CRAN.R-project.org/package=DHARMa

Rodríguez-Sánchez, F. (2023). DHARMa.helpers: Helper functions to check models not (yet) directly supported by DHARMa.